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Abstract 

In this work, we compare three turbulence models used to parameterize the oceanic 
boundary layer. These three models depend on the bulk Richardson number, which is co- 
herent with the studied region, the West Pacific Warm Pool, because of the large mean 
shear associated with the equatorial undercurrent. One of these models, called R224, is new 
and the others are Pacanowski and Philander's model (R213 model) and Gent's model (R23 
model). The numerical implementation is based on a non-conservative numerical scheme. 
The following (three criteria) are used to compare the models: the surface current intensity, 
the pycnocline's form and the mixed layer depth. We initialize the code with realistic velocity 
and density profiles thanks the TOGA-TAO array (McPhaden, 1995, [21]). In case of static 
instability zone on the initial density profile, only the R224 model gives realistic results. Af- 
terwards, we study a mixed layer induced by the wind stress. In this case, the R224 results 
and the Pacanowski and Philander's results are similar. Furthermore, we simulate a long 
time case. We obtain a linear solution for all models that is in agreement with Bennis and 
al p. 

Summary 0.1 Keywords: vertical mixing, Richardson number, mixed layer. 



1 Introduction 



The presence of an homogeneous layer near the surface of the ocean has been observed since 
a long time. This layer presents almost constant profiles of temperature and salinity. We 
distinguish the mixing layer from the mixed layer (Brainerd and Gregg, 1995, [3]). The mixing 
layer is actively mixed by surface fluxes. It is a depth zone where the turbulence is strong. 
The mixed layer is a maximum depth zone which has been mixed in the recent past by surface 
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fluxes. This layer includes the mixing layer. In this work, we study the mixed layer. The bottom 
of the mixed layer corresponds either to the top of the thermocline, zone of large gradients of 
temperature, or to the top of the zone where haline stratification is observed (Vialard and 
Delecluse, 1998, [30]). Some attempts to describe this phenomenon can be found for example 
in Defant [5] or in Lewandowski [18]. The effect of the wind-stress acting on the sea-surface is 
considered to be the main forcing of this boundary layer. Observations in situ were completed 
by laboratory experiments (Deardoff, 1969, [I]) and more recently by numerical modelizations 
of the mixed layer. A historical of these different approaches can be found in Kase ([14]). Kase 
quoted the model of Kraus and Turner (1967, [T5]) as the first applied mixed layer model which 
is a bulk one. There are also first order models as the Pacanowski and Philander model (called 
PP model, 1981, [26]) and the Large and Gent model (called KPP model, 1994, [E]). The 
second order models have been developed by Mellor and Yamada (called MY model, 1982, |22j ) 
and Gaspar and al (1990, [7])- In this work, we focus our attention on the first order models. 
According to observations, the thickness of the mixed layer can vary between ten meters and a 
few hundred of meters, depending on the latitude (Boyer de Montegut and al, 2004, [53]). The 
mixed layer depth is difficult to determine. Traditionally, there are two types of criteria for this 
determination: the density difference criterion which will be used in this work and the density 
gradient criterion. The first one consists in estimating the mixed layer depth whose the density 
is about the surface density increased by 0.01 kg.m~ 3 '. This criterion is often used in the tropical 
area as in Peters and al, 1988 [27] . So, the mixed layer is the one where the density is inferior to 
the surface density to 0.01 kg.m~ s . The second criterion is based on specifying density gradient. 
The mixed layer depth is where the density gradient is equal to the specified value. 

Mixing processes are intense in the homogeneous boundary layer, much weaker near the py- 
cnocline or in the deep ocean. The effects of local fluxes of heat and momentum across the 
sea-surface induce turbulence and then mixing processes in the surface layer. The dynamical 
response of the ocean, and especially small scale turbulence, has a significant effect on mixing 
processes. This is true particularly in tropical areas. The rate of kinetic energy dissipation and 
the vertical turbulent fluxes of heat and mass cannot be measured directly, but they can be 
deduced from measurements of vertical temperature gradients and horizontal velocity [29]. Such 
experiments shown that turbulent dissipation is higher near the equator than in low latitudes 
|10j . [llj . [BJ. Therefore in order to modelize the mixed boundary layer and to represent cor- 
rectly the mixing processes, it is necessary to define a parametrization of turbulent diffusion. The 
mixed layer being strongly dominated by vertical fluxes, attention is drawn on vertical mixing 
which requires a closure model in order to represent the Reynolds stress. Recently, Goosse and 
al [9] studied the sensitivity of a global model to different parametrizations of vertical mixing. 
They insist on the crucial role of mixing in the upper oceanic layers: it has a direct impact on 
the sea surface temperature and then on ice evolution, but affects also the vertical profile of 
velocity by redistributing the momentum [2]. 

Classically, mixing parametrizations consist in the definition of turbulent eddy viscosity and 
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diffusivity coefficients. These coefficients can be chosen either as constants or as fixed profiles 
|13j . This is a simple but crude parametrization since variations of mixing with time and 
location are forbidden. A more suitable method is to define these coefficients as functions of 
processes governing the mixing. In tropical oceans vertical stratification and velocity shear 
are natural parameters since, following Philander [29], one of the reasons for higher turbulent 
dissipation near the equator than in low latitudes, is the large vertical shear observed in tropical 
currents. Pacanowski and Philander [2B] proposed a formulation for eddy viscosity and diffusivity 
coefficients depending on the Richardson number which represents the ratio between buoyancy 
effects and vertical shear. The Richardson number dependent formulations allow strong mixing 
in high shear regions with low stratification, and low mixing elsewhere. It has to be noted that 
stratification tends to reduce turbulence and therefore mixing processes. 

In this paper, we study first order models. We focus on the behavior of Richardson number de- 
pending models: the basic one proposed by Pacanowski-Philander (PP model) and two variants 
used in Gent [Sj, Blanke and Delecluse [2J,Goosse and al [9]. In this study, the PP model is 
called R213 model and the Gent model is called R23 model. The last model is named R224. The 
studied region is the West Pacific Warm Pool. At this location, the modelization based on the 
Richardson number is realistic because of the large mean shear associated with the equatorial 
undercurrent. Peters and al [27] have shown that the PP model underestimates the turbulent 
mixing at low Richardson number, while overestimating the turbulent mixing at high Richardson 
number in comparison with turbulent measurements. They found that the simulated thermo- 
cline is much too diffused compared with observations. Furthermore, this scheme overestimates 
the surface current intensity. Also, the simulated equatorial undercurrent is too shallow com- 
pared to observations (|12].[21E].|19|). Halpern and al [12] compared the PP scheme and the MY 
scheme. At the equator, the current and temperature simulated by the PP scheme are more 
realistic than those simulated by the Mellor and Yamada scheme. The Gent model, one of three 
studied ones, gives realistic results in the West Pacific Warm Pool. This model simulates a sharp 
thermocline which is in agreement with the observations [B]. Furthermore, it gives good results 
for the annual average SST (Sea Surface Temperature) at the equator. 

In this paper, we investigate in case of each model, the simulated mixed layer depth, the form 
of the simulated pycnocline (sharp or diffuse) and the intensity of the simulated surface current. 
The numerical implementation is based on a non-conservative numerical scheme. We initialize 
the code with realistic velocity and density profiles. These profiles come from to the TOGA- 
TAO array (McPhaden, 1995, [21]). The geographic location is 0°N, 165°-E which is found in 
the West Pacific Warm Pool. In the first part, we study a mixed layer induced by the wind 
stress. Then, we initialize the code with a density profile showing a static instability zone. At 
last, we simulate a long time case. 
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2 Modelization of the mixed layer 
2.1 Setting of the problem 

Variables used in this paper to describe the mixed layer are statistical means of the horizontal 
velocity and of the density denoted by (u, v, p) . In the ocean, the density is a function of 
temperature and the salinity through a state equation. So, we consider the density as an idealized 
thermodynamic variable which is intended to represent temperature and salinity variations. 
The formation of the mixing layer is a response to sea-air interactions: wind-stress, solar 
heating, precipitations or evaporation acting on the sea surface. The variability of temperature 
is considered as essential to understand the response of the ocean. For example, the existence 
of a sharp pycnocline is a well known feature of the tropical areas. The thermal inertia of 
the water column is linked with the depth of the pycnocline, which influences the sea-surface 
temperature. The role of the haline stratification, sometimes considered as less important, has 
been recently evidenced by Vialard and Delecluse [30J: a numerical modelization of the tropical 
Pacific produces a "barrier layer" depending on surface forcings and large scale circulation. The 
term "barrier layer" refers to the water column located between the bottom of the mixed layer 
and the top of the pycnocline. It is present when the isohaline layer is shallower than the 
isothermal layer and then the depth of the mixed layer is controlled by the salinity. 
The model studied hereafter is not expected to describe all the phenomena occuring in the mixed 
layer. Its purpose is only a better understanding of a classical closure model. Therefore we shall 
use simplified equations governing the variables u (zonal velocity) ,v (meridian velocity) and p 
(density). 

The mixed layer being strongly dominated by vertical fluxes, we shall suppose that u,v and p 
are horizontally homogeneous and we so obtain a one-dimensional modelization. The Coriolis 
force will be neglected, which is valid in tropical oceans. Therefore, the equations governing the 
mixing layer are 



where u',v',w',p' represent the fluctuations of the horizontal velocity, vertical velocity and the 
density. The notation ( ) signifies that the quantity is statistically averaged. Equations (1) are 
the classical equations corresponding to a modelization of a water column. Such equations can 
be found in [T7] as the equations of the boundary layer. Equations (1) are not closed and then 
the vertical fluxes appearing in the right-hand side have to be modelized. 

We study in this paper the behavior of a very classical closure modelization that uses the concept 
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of eddy coefficients in order to represent turbulent fluxes. So we set 

du 



[u w') 



(v'w') 



{p'w') 



V2 



dz ' 
dv 

dz' 
dp 

dz' 



Coefficients v\ and v 2 are called vertical eddy viscosity and diffusivity coefficients and will be 
expressed as functions of the Richardson number R defined as 



R 



dp 

dz 



(H) 2 + (ff) 2 ' 



where g is the gravitational acceleration and p a reference density (for example p = 1025 kg.m ). 
The set of equations, initial and boundary conditions governing the mixing layer can now be 
written 



du d ( du 



dt dz 
dv d 



dt dz V 1 dz 



dz 
dv 



dp d 



dt dz \ U2 dz 



dp 



= 0, 
0. 
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u = Ub, v = Vb, p = Pb at the depth z 



-h. 



du p aTr dv p aTr dp 
v\— = — V x , v\— = — V y , v-i-^- = Q at the surface z = U, 
dz p dz p dz 

u = uq, v = vq, p = Pq at initial time t = 0. 



In system (2), the constant h denotes the thickness of the studied layer that must contain the 
mixing layer. Therefore the circulation for z < —h, under the boundary layer, is supposed to be 
known, either by observations or by a deep circulation numerical model. This justifies the choice 
of Dirichlet boundary conditions at z = —h, Ub, Vb and p h being the values of horizontal velocity 
and density in the layer located below the mixed layer. The air-sea interactions are represented by 
the fluxes at the sea-surface: V x and V y are respectively the forcing exerced by the zonal wind- 
stress and the meridional wind-stress and Q represents the thermodynamical fluxes, heating 

I 1 2 I 1 2 

or cooling, precipitations or evaporation. We have V x = Cd \u a \ and V y = Cd \v a \ , where 
U a = (u a ,v a ) is the air velocity and Cx>(= 1,2.10 -3 ) a friction coefficient. 



We study hereafter three different formulations for the eddy coefficients v\ 
fo(R)- Functions f\ and fa can be defined as 



fi (R) and u 2 
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Formulation (3) corresponds to the modelization of the vertical mixing proposed by Pacanowski 
and Philander [26]. They proposed for coefficients a^/^and a 2 the following values: a\ = 
1.10 -4 , 0i = 1.10 -2 , a 2 — l-10 _5 (units: m 2 s~ l ). This formulation has been used in the OPA 
code developed by Paris 6 University [2], [20] with coefficients a\ = 1.10" 6 , f3 1 = 1.10~ 2 , a 2 = 
1.10~ 7 (units: m 2 s _1 ).The selection criterion for the coefficients appearing in these formulas was 
the best agreement of numerical results with observations carried out in different tropical areas. 
A variant of formulation (3), proposed by Gent [8J, is 



— 2, f2(R) = a 2 + 



f 1 (R) = a 1 + ^ 7 „ , f 2 (R) = a 2 + - (4) 



with a x = 1.10~ 4 , /?! = 1.10 -1 , a 2 = 1.10" 5 , @ 2 = 1-10 -1 (units: m 2 ^" 1 ). A formulation 
similar to (4) when replacing 10-R by 5R and varying the values of the coefficients a±, a 2 between 
the surface and the depth 50m is used in [9]. 

In this paper, we will also study the properties of another formulation close to formula (3): 



h h{R)=a2+ tm = a2+ * + h 

(l + 5i?) 2 ' h [ ) a2+ (l + 5i?) 2 ° 2 + (1 + 5.R) 2 + (1 + BR) 4 



with ax = 1.10" 4 , 0x = 1.10~ 2 , a 2 = 1.10" 5 , 2 = 1.10" 3 (units: m 2 S - 1 ) . 

Eddy viscosity vx defined by (5) is the same as the coefficient given by Pacanowski and Philander. 

The definition of the eddy diffusivity coefficient v 2 differs by the exponent of the term (1 + 5R). 

In formulas (3) to (5), the eddy coefficients vx and v 2 are defined as functions of the Richardson 
number R through the terms (1 + 7-R) n appearing at the denominator. Hereafter, these three 
formulations will be denoted respectively by R213, R23 and R224 where R signifies Richardson 
number and the integer values are the exponents of (1 + "yR) in the definitions of vx and v 2 . 
Eddy coefficients defined by relations (3) to (5) all present a singularity for a negative value 
of the Richardson number R = —0.2 or R = —0.1. We have plotted in Figure la the curves 
vx = fi(R)- In formulations (3) or (4) the coefficient of eddy diffusivity u 2 becomes negative for 
values of R lower than —0.2 or —0.1, and therefore the models are no more valid. The curves 
v 2 — fi (R) obtained with formulations (3) (4) (5) and are plotted in Figure lb. 
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Figure la: Viscosity {y\ = fi(R)) for all models 
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Figure lb: Diffusivity {y<i = f2{R)) fo r all models 

Problem (2) coupled with one of the definitions (3) to (5) for eddy coefficients retains the vertical 
shear and buoyancy effect which are two important processes for the generation of the mixed 
boundary layer, especially in tropical areas. 
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3 Numerical Experiments 



3.1 Finite Difference Scheme 

We want to resolve numerically the system ([2]). We replace the continuous variables (u, v, p, v\, zv 2 ) 
by discrete variables (u™, vf,p™, (z^i)™ , (^2)") which are the approximate solutions at time nAt 
(with n = 1,2, ...,N) and at points (i — NI)Az (with i = 1,2, NI). We discretize the 
ID-domain in z-levels where z is the vertical coordinate. 

We use a second-order central difference scheme for the second space derivative and a first-order 
backward difference scheme for the first space derivative. These previous schemes can be written 
as: 



Second-order central difference scheme: ( 



fd 2 u\ n+1 _ u#l - 2< +1 + <+/ 



2 



V dz 2 J i Az 

/du\ n+1 ^ +1 -»' n+1 
First-order backward difference scheme: ( — — 



\dz J i Az 

The grid spacing, Az, is equal to 5 m in cases 2 and 3 and equal to 1 m in case 1. The time step, 
At, is equal to 60 s. In time, we use implicit velocities and implicit density. The viscosity {u\) 
and diffusivity (z/2) are explicit. The basin depth is 100 m. The boundary conditions are treated 
with a first-order backward difference scheme. We choose Neumann boundary conditions at the 
surface and Dirichlet boundary conditions in z = —h. 

The numerical scheme is the following: 
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and the Neumann boundary conditions at the surface are computed as: 

/„.n+l „,n+l\ / „,n+l „,n+l\ / n n+l n n+l\ 

^■O*-^)-^ 

Moreover, the residual values are computed as: 

/ NI NI NI \ 1 I 2 

r n = E K +1 - «?l a + E K +1 - <\ 2 + E \p T i +l - ^l 2 • ( ? ) 

\i=l i=l i=l / 



(6) 
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Our numerical scheme is non-conservative but it produces the same results as those given by 
a conservative one. Furthermore, the residual values decrease monotonically to zero as the 
flow reaches the steady states, better than the conservative scheme where the behaviour of the 
residual is more erratic. 

3.2 Numerical Results 

In this section, we aim to study an Equatorial Pacific region called the West-Pacific Warm Pool. 
So, we initialize the code with data from the TAO array (McPhaden [21]). It is located at the 
equator between 120° E and lSO ^. The sea temperature is high and almost constant along the 
year (28 — 30° C). The precipitation are intense and hence the salinity is low. In the first part, 
we study a mixed layer induced by the wind stress. Then, we consider a case where the code 
is initialized by a density profile showing a static instability zone. Finally, we simulate a long 
time case. All our numerical results are grid-size (Az) and time step (At) independent, in the 
sense that they remain practically unchanged as Az and At decrease. 

3.2.1 The initial data 

We use data available from the Tropical Atmosphere Ocean (TAO) array (McPhaden [SJJ ). The 
TAO project aims to study the exchange between the tropical oceans and the atmosphere. The 
TAO data have being very used in numerical simulations. The velocity data comes from the 
ACDP (Acoustic Doppler Current Profiler) measurements. To obtain the appropriate profiles, we 
interpolate the data by a one-order linear interpolation. We initialize our code by these profiles at 
0°iV, 165° E and we obtain the results below. The buoyancy flux is equal to — 1.10~ 6 kg.m~ 2 .s" 1 
(~ — 11 W/m 2 ) in all cases. This heat flux is in agreement with Gent [8] because the heat flux 
must be in the range [0 W/m 2 ; 20 W/m 2 } between 140° E - 180° E and 10° N - 10° S. 

3.2.2 Case 1: A mixed layer induced by the wind stress 

We use, in this case, as initial data, the TAO's data for the time period between the 15th June 
1991 and the 15th July 1991. The initial profiles are displayed on figure 2. The initial zonal 
velocity profile presents a westward current at the surface and, below it, an eastward under- 
current whose maximum is located about 55 m. Deepest, we observe a westward undercurrent. 
The initial density profile does not display a mixed layer. 
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Figure 2: Initial zonal velocity, meridian velocity and density profiles (from left to right). 

The buoyancy flux is equal to — 1.10 -6 kg.m~ 2 (~ —HW/m 2 ). The model is integrated 
for 48 hours. The grid spacing is equal to 1 meter and the time step is equal to 60 seconds. 
The chosen wind stress is stronger than in reality because we want to simulate a mixed layer 
induced by the wind stress. These values correspond to an another period in the studied year. 
The zonal wind is equal to 8.1 m/s (eastward wind) and the meridional wind is equal to 2.1 m/s 
(northward wind). 

The numerical results are displayed on figure 4. The final density profile displays a similar mixed 
layer for R213, R224, R23 models. Furthermore, the pycnocline simulated by R224, R213 and 
R23 models are similar. 

The deep flow (60 — 100 m) are the same for all models because the surface fluxes are not strong 
enough to affect the deep water. 

The R213 and R224 surface current are of the same kind. The R23 model underestimates this 
current. We notice an increase in zonal and meridian surface current in comparison with the 
initial profiles which concord with a south-westerly wind. The surface current behavior can be 
explained by the viscosity and diffusivity values. The final diffusivity and viscosity are displayed 
on figure 3 for all models. We observe the order described below. 

("1)23 > ("1)213 > ("1)224 

(1/2)23 > ("2)213 > ("2)224 

The R23 model has a strongest viscosity and diffusivity. Therefore, the R23 surface current is 
lower than the other surface currents. As the R224 and R213 viscosities and R224 and R213 
diffusivities are of the same kind, the R224 and R213 surface currents are similar. However, the 
R224 surface current is slightly stronger than the R213 surface current. 
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Figure 3: Final viscosity (left-hand side) and final diffusivity (right-hand side) for all models. 
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Figure 4: Comparison of three turbulence models. We can see, on the top row, the entire 
water column and on the bottom row, the shallow flow. We have plotted, from left to right, the 
finals zonal velocity profiles, the final meridian velocity profiles and the final density profiles. 
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Figure 5: Comparison of residual values. 
On figure 5, we notice a monotonic numerical convergence to the steady state for all models. 

3.2.3 Case 2: A static instability zone in the initial density profile 

The code is initialized with the 17th November 1991 data. The initial zonal velocity profile 
displays an eastward current whose maximum is located about 55 m. The initial meridian 
velocity profile displays a southward current whose maximum is located about 20 m. The initial 
density profile displays a static instability zone between —30 m and —50 m. Notice that there 
is a similar instability zone for the following days: 16th November, the 19th November, the 
20th November, the 21th November and the 22th November 1991. Here, we observe a seventy 
meters deep mixed layer. However, this mixed layer is not homogeneous. The initial profiles are 
displayed on figure 6. 




Figure 6: The initial profiles for zonal velocity, meridian velocity and density (from left to 

right). 

On figure 7, the initial Richardson number, called BP, is near to —0.2 for z = —45 m. The 
R213 and R224 models have infinite viscosity and diffusivity for R = —0.2 (see Figures la and 
lb). Hence, the R213 diffusivity (see figure 9) and the R224 diffusivity (see figure 8) are large 
for this depth. The initial Richardson number (see figure 7) is inferior to —0.2 for z = —35 m 



12 



and z = —50 m. Therefore, the R213 diffusivity is negative for this depths (see Figure lb). 
In the range [—35 m, —50 m], the initial richardson number is inferior to —0.1 and hence the 
R23 diffusivity is negative (see Figure lb). On figure 9, the negative diffusivities are marked 
by a point for R23 and R213 models. The R224 diffusivity (see figure 8) is always positive. In 
physical reality, negative diffusivity does not exist. The diffusivity was estimated by Osborn 
and Cox [25] with measurements of very small scale vertical structure. They have shown that 
in our studied region, the diffusivity was in the range [1.10 -2 cm 2 .s _1 , 1.10 3 cm .s ]. So, we 
can not use R213 and R23 models in this case. 
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Initial Richardson number Initial Richardson number 

Figure 7: The initial richardson number for different depths. We can see a entire water 
column on the left-hand side and the flow between -30m and -55m on the right-hand side. 
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Figure 8: Initial diffusivity for formulation R224. 
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Figure 9: Initial diffusivity for formulations R213 (left-hand side) and R23 (right-hand side). 
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Figure 10: Zoom on the R213 negative values marked by an asterisk. 



The model is integrated for 48 hours. The grid spacing is equal to 5 meters and the time step 
is equal to 60 seconds. The zonal wind is equal to 11.7m.s _1 (eastward wind). The meridional 
wind is equal to 0.4 m.s~ l (northward wind). The buoyancy flux is equal to — 1.10" 6 kg.m~ 2 
(~ -llW/m 2 ). 

The results are displayed on figure 11. The R224 model produces a homogeneous seventy meters 
deep mixed layer. This layer is homogeneous because we have applied a negative buoyancy flux 
at the surface which stabilizes the flow. 

We notice an increase in the zonal surface current, in comparison with the initial surface current, 
which is in agreement with an eastward wind at the surface. The northward wind at the surface 
is too weak to really modify the initial meridian surface current. 
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Figure 11: Zonal velocity profile (left position), meridional velocity profile (medium position) 
and density profile (right position) in case of the R224 formulation. 
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Figure 12: Evolution of residual values in case of the R224 formulation. 

On figure 12, the residual values display a good numerical convergence for R224 model (solid 
line) . 



3.2.4 Case 3: A long time case 



We simulate a long time case. So, the model is integrated for 10000 hours. The initial profiles rep- 
resent the ocean mean state on June 17, 1991. The buoyancy flux is equal to — 1.10" 6 kg.m~ 2 
(~ —11 W/m 2 ). The surface zonal wind is equal to 5.4 m.s -1 (eastward wind) and the surface 
meridional wind is equal to 0.9 m.s" 1 (northward wind). 

The initial zonal velocity profile (see Figure 13) displays two eastward currents located at the 
surface and around to — 70 m as well as two westward currents located around to —45 m and 
—90 m. The initial meridian velocity profile (see Figure 13) displays several southward currents 
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and northward currents. The main southward current is located around to —55 m while the 
main northward current are located around to —90 m and at the surface. Furthermore, on 
initial density profile (see Figure 13), we notice a thirty five meters deep mixed layer according 
to Peters and al's [281 criterion. 
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Figure 13: The initial profile of zonal velocity, meridional velocity, density (left to right). 

The numerical results are displayed on figure 14. The three turbulence models give a linear 
profile for the simulated time. This fact corroborates the existence of a linear equilibrium 
solution obtained by Bennis and al pQ. 
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Figure 14: Comparison of zonal velocity profiles (left position), meridional velocity profiles 
(medium position) and density profiles (right position). 

The residual values are displayed on figure 15. Notice that the numerical convergence to steady 
states is good. 
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Figure 15: Comparison of residual values. 

As the equilibrium solutions (zonal velocity u e , meridian velocity v e and density p e ) are unique 
in case of the R224 model (see Bennis and al, 2007 [I]), we can compare the theoretical solutions 
with the numerical solutions. We remember hereafter the theoretical equilibrium solutions which 
satisfy the following system: 

I (* <*•>£)=»■ I (* 

The theoretical steady-state profiles for velocities and density are obtained by integrating the 
previous system with respect to z, taking account the boundary conditions at z = —h: 

u e (z) = u b + Vxpa (z + h), v e (z) = v b + V f a (z + h), f( z )=p b+ -Q-(z + h). 
Poh \R ) Po/i ( R ) h \R ) 

The equilibrium richardson number {R e ) can be interpreted as the intersection of the curves 

pl(v x 2 + v v 2 ) 

-.In case of R224 model, the graph 



k(R) 



(a my 



and h (R) = CR with C 



f 2 (R) - gQp 

of function k and h for Q < is plotted on figure 16. 
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Figure 16: Equilibrium richardson number - R224 model. 
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So, we compute the theoretical solutions with the equilibrium richardson number (Re = 0.063935) 
given by the equation h(R) = k(R). The theoretical and numerical equilibrium solutions are 
displayed on figure 17. 
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Figure 17: Theoretical and numerical equilibrium solutions - R224 model. 
We obtain similar profiles and the small differences are due to numerical scheme which we use. 

4 Conclusion 

4.1 Code validation 

We observe, in case 1, a mixed layer induced by the wind stress that is a physical phenomenon. 
Moreover, we notice, in all cases, an increase in the zonal surface current when we apply an 
eastward wind (u > 0) at the surface. In the same order, an efficiency northward wind (v > 0) 
at the surface causes an increase in the meridional surface current. These results are in agreement 
with the physical reality. Furthermore, the residual values correctly converge to steady states, 
corresponding to rather high levels of turbulent diffusion. Theses previous observations validate 
our code that gives a good representation of the physical reality. 

4.2 Comparison of three turbulence models 

Our comparison is based on three criteria: the mixed layer depth, the surface current intensity 
and the pycnocline's form. The mixed layer depth is obtained with a density difference criterion. 
This difference is equal to 0.01 kg.m~ 3 . We use the surface current values to determine the 
surface current intensity. The density gradient is used to determine the pycnocline's form. We 
summarize the results hereafter. 

- In the particular case 1, the R213, R23 and R224 models give a same mixed layer depth. In 
terms of surface current, those simulated by R224 and R213 models are of the same kind. 
However the R224 surface current is slightly stronger than the R213 surface current. The 



18 



R23 model underestimates this current. Furthermore, the pycnocline's form are similar 
for R213, R23 and R224 models. 

- In the particular case 2, the R213 and R23 models have negative diffusivities at the initial 
time. Therefore, these models are not physically valid in this case. This problem comes 
from static instability zone in the initial density profile. Hence, we can not use these 
models in this case. Then, only R224 is valid. The R224 model produces a homogeneous 
mixed layer. So, we do not compare the pycnocline's form since we have only one mixed 
layer. 

- Furthermore, the long time profiles are linear for all models that is in agreement with 
Bennis and al PQ. In case of R224 model, the theoretical and the numerical profiles are 
similar. 

Globally, the R224 model has the same behavior as the Pacanowski and Philander model (R213 
model) and we can use it in more situation. For example, R224 model can be used in almost 
all of convective cases. Moreover, the R224 numerical equilibrium solution is in agreement with 
the R224 theoretical equilibrium solution. 
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